Noname manuscript No. 

(will be inserted by the editor) 



Transition to chaos and escape phenomenon in two degrees 
of freedom oscillator with a kinematic excitation 

Marek Borowiec and Grzegorz Litak^ 



Received: date / Accepted: date 



Abstract We study the dynamics of a two-degrees-of-freedom (two DOF) non- 
linear oscillator representing a quarter-car model excited by a road roughness 
profile. Modelling the road profile by means of a harmonic function we derive 
the Melnikov criterion for a system transition to chaos or escape. The analyti- 
cally obtained estimations are confirmed by numerical simulations. To analyze the 
transient vibrations we used recurrences. 
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1 Introduction 

The dynamics of a quarter-car model is governed by the road profile excitations 
and nonlinear suspension characteristics [UGHS] • These elements were discussed 
separately using simplified models and jointly by considering the more realistic de- 
scription of a vehicle motion. For instance Verros et al [1] proposed a quarter-car 
model with a piecewise linear dynamical characteristics. According to the adapted 
control strategy, the damping coefficient switched between two different values. 
Gobbi and Mastinu [2] considered a two DOF model to derive a number of ana- 
lytical formulae describing the dynamic behaviour of passively suspended vehicles 
running on randomly profiled roads. Their linear model approach was generalized 
by Von Wagner [3] who determined a corresponding high-dimensional probability 
density by solving the Fokker-Planck equations. Finally, in papers [HEM] dynam- 
ics, bifurcations and appearance of chaotic solutions were discussed. 

However the main issues of vehicle dynamics studies were unwanted and harm- 
ful vibration responses generated by vehicle as an effect of a rough surface road 
profile kinematic excitation p1HlF7l[8l[9l ll0llllj Thus the efficient reduction of them 
is still a subject of research among automotive manufacturers and research groups 
[14U15] . Turkay and Akcay [7] considered constraints on the transfer functions from 

Department of Applied Mechanics, Lublin University of Technology, Nadbystrzycka 36, PL- 
20-618 Lublin, Poland, 

# corresponding author: tel.: +4881 5384573 
fax: +4881 5384205 E-mail: g.litak@pollub.pl 



2 



Marek Borowiec and Grzegorz Litak* 



the road disturbance to the vertical acceleration, the suspension travel, and the 
tire deflection are derived for a quarter-car active suspension system using the 
vertical acceleration and/or the suspension travel measurements for feedback. The 
recent experimental and theoretical studies involved many new applications of ac- 
tive and semi-active control procedures 12, 13, 141 1151116] . Consequently, previous 
mechanical quarter car models [3l ll0l[TT] ll5] have been re-examined in the con- 
text of active damper applications. Dampers based on a magnetorheological fluid 
with typical hysteretic characteristics have a significant expectation for effective 
vibration damping in many applications [4l[T7l[T8l[T9l[20ll21U22U23] . 

Recent efforts have been also focused on studies of the excitation of an auto- 
mobile by a road surface profile with harmful noise components [31 1241125] . Noise 
like chaotic vibrations, appearing due to the system nonlinearities, have been in- 
vestigated in simplified single DOF models [4l|T9"||6] • These papers follows the rich 
literature on an escape phenomenon in the symmetric and nonsymmetric Duffing 
or Helmholtz potentials [26,27,28,29,30_, where critical system parameters were 
determined. Note that the condition for escape could be applied as a criterion of 
fractality in basins of attraction and also as a transition to chaos [?] . On the other 
hand, more sophisticated models of vehicle dynamics, namely half-car and full-car 
models, in the context of nonlinear response including chaotic solutions have been 
studied by Zhu and Ishitobi and also by Wang et al. 31,32,33 . 

To study a transition to the chaotic region and the corresponding critical pa- 
rameters in low dimension dynamical system analytically, the Melnikov theory 
34,35,36,37,38 is often advocated. The application of this approach to a simple 
quarter-car model has been recently proposed by Li et al., and Litak et al. [H 
119] , In the above papers, a single DOF model was used because of its simplicity. 
Consequently, the analytic consideration included also multiple scales analysis and 
harmonic balance [201139] . The present paper is the continuation of the previous 
studies with an extension to a more realistic two DOF model, which includes the 
sprung and unsprung masses (Fig. 1). 
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Fig. 1 Two-degree-of-froedom quarter-car model (y = x± — x^)- 



Our article is organized in 4 sections. After introduction in the present section 
(Sec. 1) we present the model and discuss possibility of a global homoclinic bifur- 
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cation (Sec. 2). By the reduction of dimension we introduce the foundations of the 
Melnikov approach to the model. In this section (Sec. 2) we obtain the principal 
result as a critical curve which define the system parameters regions of regular 
and non-periodic (chaos, transient chaos or escape) behaviour. The simulation re- 
sults, illustrating that transition results are also shown. Further results including 
recurrence analysis confirming the theoretical predictions are provided in Sec. 3. 
Finally in Sec. 4 we end up with conclusions and final remarks. 



2 The model and global bifurcations 

We start the analysis from the two DOF model presented in Fig. 1. The dynamics 
of vehicle excited by the road profile quarter car model is governed by nonlinear 
suspension characteristics. To examine a transition to the chaotic regime of vibra- 
tions, the Melnikov theory has been recently proposed [4"lll9|. In this perturbation 
approach the authors of previous works used a single DOF model. 

In the present note we go beyond this assumption by considering an extension 
of a vehicle model with the defined unsprung and sprung masses (Fig. [1^) . The 
differential equations of motion for both masses have the following form: 

d 2 ci , d d . ai , , , , 

-T7o x i H (^7 Xl ~ TT 12 ) H (X1-X2) (1) 

dt- mi at at mi 

+ ^{x 1 -x 2 f = 
mi 

d 2 c 2 , d d . ci . d d . , . 

-77^2 H (—X2--T-X0) (-77 XI - — X2) (2) 

at- m,2 at at m,2 at at 

H {X2-X0J (-31-352) (351-352) =0 

m2 m2 m,2 

where mi, m2 denote the corresponding sprung and unsprung masses (Fig. [I} 
xq = acosujt and describes the harmonic corrugation of a road profile, c,, ai (for 
i = 1, 2) & /3i are damping and stiffness coefficients, respectively. 

Now we define the new variable y = x\ — X2 of the relative motion: 

■■ Cl Ql ft 3 .. ,„s 

y + t — y-\ y-\ y =x 2 . (3) 

mi mi mi 

To adjust the above equations to application of the higher dimensional Mel- 
nikov approach [40] .the above equations can be approximated by introducing a 
small parameter e: 

■■ i Cl • , «1 , Pi 3 •• ( a\ 

y + t — y-\ y-\ y = -x 2 - (4) 

mi mi mi 
Consequently, the second equation (Eq. [2]) gets the new form: 

C2 . 012 OL\ Pi 3 C"2 . Oi 2 , . 

X2 H X2 H 352 = — y H y waesinwtH aecosujt. (5) 

m2 m2 m2 m,2 m2 m2 



Using Melnikov approach [3l] for y coordinate equation (Eq. 3) we study the 
Hamiltonian system of the nodal kinetic energy defined at the saddle points. In 



1 



Marek Borowiec and Grzegorz Litak* 



this way, time variability of y and y are negligible. The above equations can be 
rewritten in the dimensionless form: 



y + eCiA iy + y + A 3 y 3 = -x%. (6) 



x 2 + C 2 A 1 x 2 + A 2 Mx 2 = C!M Vmiai y + My + A 3 My 3 (7) 

mi 

— eC2Mausm{ijjt) + eA2Macos(tot), 

where Ci = ^,C 2 = A 1 = ^SjL, A 2 = ^, A 3 = ^, M = ^ are dimension- 
less parameters. For further consideration we assumed that the system parameters 
are: Ci = 0.001, C 2 = 0.5, Ai = 1, A 2 = 1, A 3 = -16, M = 5. The excitation fre- 
quency and the corresponding amplitudes used in the analysis have been fixed to 
uj = 1.5, a = 0.08 and a = 0.12. In the numerical simulations we used the sampling 
time St = 0.00418. 

To analyze a homoclinic bifurcation in the sprung mass vibration we make 
further approximation. The role of the small parameter e is to determine the het- 
eroclinic trajectory and decouple the equations of motion into separate equations 
for sprung and unsprung masses. Thus after above normalizations the equations 
can be expressed: 



v = -eCiv - y - A 3 y 3 - i 2 , (8) 
y = v, 

X2 = -C 2 x 2 - A 2 Mx 2 + M (dv + y + A 3 y 3 (9) 

— tC2ua smuit + eA2acostot) . 

Interestingly, in the limit of small e limit, X2 can be approximated as: 

x 2 = eAcos(ut + 0- ip) + x st . (10) 

Note that in the above expression x s t , generated by slowly changing terms with y 
and y, is plating a role of the static displacement while A as an amplitude. 

Thus, the unperturbed equations (for e = 0) can be obtained from the gradient 
of the unperturbed Hamiltonian H°(y,v): 

. 8H° . dH° 

^=^T' Vy = —8y-' 



where H° is defined as follows: 



H° = ^ + V(y) (12) 
The corresponding effective potential is given by the expression: 

V(y) = ^+A 3 ^. (13) 

Following the standard Melnikov theory [31M51I361I571I58] . we get the hete- 
roclinic orbits connecting the two hyperbolic saddle points (coinciding with the 



maxima of the potential V(y) (Eq. 12)): y = ±\i A 
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They can be expressed analytically as: 



A 3 1 tanh 



(*-*o) ; 



dy*(t) 
dt 



V2 



2 cosh 2 (t-to) 



>2A: 



(14) 
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Fig. 2 Potential of the restore force - V(y) (Eq. 12) for different A3 (a) and the corresponding 
heteroclinic orbits (b). 



After perturbations the heteroclinic orbits the stable and unstable manifolds 
are calculated. Because of perturbations they are detracted (Fig. 3). As the char- 
acteristic distance between them d — > the system possibilities of mixed solutions 
(regular and escape) appear which is equivalent to irregular chaotic and chaotic 
transient solutions. Thus d = is the ideal criterion for chaos appearance. 



d 




Fig. 3 Stable and unstable manifolds of perturbed orbits terminated and started in the cor- 
responding saddle points, and the distance between them d (in y — v plane). 



Formally, the distance d between perturbed stable and unstable manifolds are 
proportional to the Melnikov integral (Fig. 3) d ~ M(ta) [34,35,36 which can be 
written: 



/+00 
h(y*(t),v*(t))Ag(y*(t),v*(t))dt (15) 
-00 

where A denotes a wedge product, the differential form h is the gradient of the 
unperturbed Hamiltonian and g is related to the perturbation part, and to is an 
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1.5 2 2.5 3 3.5 4 4.5 , w 1.0 1.5 2.0 

Fig. 4 The Melnikov criterion, a = r;Ci versus a; for system parameters where A3 = —16, 
M = 5, A2 = 1.5, Ci = 0.001 C2 = 0.5. The curve separate the region of regular solutions 
(below the curve) from the chaotic and escape solutions (above the curve). Note that Fig. 4b 
magnifies the marked region in Fig. 4a, where the two indicated points symbolizes parameters 
used for numerical calculation points above and below the critical curve for a = 0.08 at the 
point T', and a = 0.12115 at the point '2' ( w = 1.5 for both cases). 
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Fig. 5 Time series (displacements of the sprung mass) of regular (a) a and irregular (or chaotic 
transient) (b) solutions for the same system parameters as in points '1' and '2' in Fig. 4a and 
b, respectively (the sampling time St = 0.00418). Note, the larger amplitude and additional 
modulation in (b). 



integration constant. Both forms are defined on unperturbed heteroclinic orbit 
stable and unstable manifolds W st{unst) = (y* s t {uns t),Kt(,unst))- 

h = vdv + (-y - A 3 y 3 )dy (16) 
g = {Alo 2 cos(ajt + <f> - 9) - C\v)Ay. 

Thus, after substitution the above forms to the Melnikov function M(to) (Eq. 
14) we get: 

M(t )=M f (-dv* +Aui 2 co8(u)t + <t>-9)\ v*dt (17) 

J — oo 

The condition for a global homoclinic transition, corresponding to a horse-shoe 
type of stable and unstable manifolds possible cross-section, can be written as: 
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V-M(to)=0 and ^ . (18) 



dt 

Consequently, the critical parameter r\ = ajC\ is now 



I {MM - + CW. sinh , _^ , (19) 



The results of the above analysis are presented in Fig. 4a. The curve separate 
the region of regular solutions from the chaotic and escape ones. The black points 
represents the parameters used for numerical simulations. To show this region 
of parameters Fig. 4b magnifies the surrounding area. The corresponding time 
histories are presented in Fig. 5a and b. Note that in purpose of simulations we 
used Eqs. 7-8 with e = 1. 
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Fig. 6 Phase portraits (displacements versus velocity of the sprung mass) of regular (a) and 
irregular (b) solutions for the same system parameters as point T' and '2' in Fig. 4a and b, 
respectively 



Note that Fig. 5a shows the mono frequency while Fig. 5b more complex re- 
sponses ( mult i- frequency or irregular). This is clearly visible in Fig. 6a and b where 
we show the corresponding phase diagrams. In Fig. 6a one can see a single line 
while in Fig. 6b lines are split into characteristic three lines patterns. It is worth 
to mention that the Melnikov criterion (Eq. 18 and Fig. 4) specifies the global 
transition associated with the destruction of borders between basins of attractions 
belonging to different solutions. In our case one of basin is related to escape from 
the potential well. The irregular solution denoted as no. 2 in Fig. 4, (see also Figs. 
5b, 6b) must be related to such an escape. Really if someone continue the calcu- 
lations for long enough time interval one observe the escape in the plot on time 
series (Fig. 7a) and on the phase portrait (Fig. 7b), respectively. 

Focusing on the above solutions we discuss the recurrence properties of numer- 
ical results with more details in the next section. 
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Fig. 7 Time series of displacement (a) and phase portrait (b) during transient vibration escape 
(corresponding time interval At = 3844.25, while sampling time St = 0.00418). 

3 Recurrence plots analysis 

The numerical solutions of the regular and transient nature can be analyzed more 
carefully by the recurrence plots. [411142] . This method is basing on the recurrences 
statistics and can be described by the matrix form R'"' 6 with corresponding and 
1 elements: 

R^' E = 0{e - |xi - Xj|) for \i - j\ > w . (20) 

where x; and xj are usually defined in the embedding space of m dimension and e is 
the threshold value. Here indices i and j denote the sampling instants. In our case 
we decided to use the m = 2 and the embedding space includes the displacements 
of sprung and unsprung masses: Xj = [y, X2] (Eqs. 7-8 with e = 1). 

Having and 1 values to be translated into the recurrence diagram as an empty 
place and a coloured dot respectively (see Fig. 9). Here w denotes the Theiler 
window used to exclude identical and neighbouring points from the analysis [H] . 
Webber and Zbilut [35] and later Marwan and collaborators [431133] developed the 
recurrence quantification analysis (RQA) for recurrence plots. 

Shortly after invention RQA were addressed to the biological and physiologic 
systems [45U46] . Recently this method has been applied for several technical sys- 
tems j47ll48il49j The first parameter of the RQA analysis defining the correlation 
function is the recurrence rate RR: 

1 N 

^=^i>r> (2i) 

i,j=l 

which calculates the number of recurrences. In our case Theiler window to = in 
order to get the consistency with correlation sum [34] . In this language e expresses 
the correlation length of the characteristic sphere radius in the embedded space. 
Note that the correlation sum is the important tool which could be used to derive 
correlation dimension D-2- For small enough e: 



D2 = logRRie) 
loge/eo 



where eq is the arbitrary length. 
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Fig. 8 The idea of states summation idea in the embedded space and the sphere of radius e. 

Table 1 Summary of recurrence quantification analysis (RQA) for m = 2 and e = 0.01 for T' 
and '2' solutions (see. Fig. 4-6). 



type of 










motion 


RR 


DET 


LAM 


Lentr 


(>i>) 


0.0348 


1.0000 


0.6078 


1.3135 


('2') 


0.0080 


0.9578 


0.4743 


1.7928 



Furthermore the RQA can be used to identify topological structures of diagonal 
and vertical lines. In its frame RQA provides us with the probability p(l) or p{v) of 
line distribution according to their lengths I or v (for diagonal and vertical lines). 
Practically they are calculated 

p( ^sf3^' (23) 

where x = I or v depending on diagonal or vertical structures in the specific recur- 
rence diagram. P £ (x) denotes the unnormalized probability for a given threshold 
value e. In this way Shannon information entropies (L^ntr) can be defined for 
diagonal line collections 

N 

^ENTR = — 

E p(0M0- (24) 

l=l m in 

Other properties, such as determinism DET and laminarity LAM: 



DET = ^'=' m " — , (25) 

£^=i^(0 
y N vp e (v) 

LAM = ^" = ""'" ^, 

where l m in arid v m i n denotes minimal values which should be chosen for a specific 
dynamical system. In our calculations we have assumed lmin — ^min — 2. 

Determinism DET is the measure of the predictability of the examined time 
series and gives the ratio of recurrent points formed in diagonals to all recurrent 
points. Note in a periodic system all points would be included in the lines. On 
the other hand laminarity LAM is a similar measure which corresponds to points 
formed in vertical lines. This measure indicates the dynamics behind sampling 
point changes. 
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Fig. 9 Recurrence plots for regular (a) and irregular (b) solutions for the same system pa- 
rameters as point '1' and '2' in Fig. 4a and b, respectively (e = 0.01). 




E 



Fig. 10 Recurrence rate RR versus threshold value e for regular '1' and irregular '2' solutions 
('1' and '2' as Fig. 4a and b, respectively). 

The results of our analysis calculated for time series '1' and '2' (see Fig. 4-6) are 
presented in Fig. 9 where we present the results of RP for e = 0.01. Note that both 
plots show regular features. However Fig. 9a is composed of full diagonals lines 
only, while in Fig. 9b each 5 line is full and between them one observes short line 
pieces or even insulated points. This would indicate multi-frequency modulated 
solution. Could be also a transient with basic regular and superimposed chaotic 
solutions. To shine this difference with more lights we show some estimated RQA 
parameters in Tab. 1. Obviously RR is smaller for '2' (Fig. 9b) as we have broken 
lines instead of full ones (Fig. 9a). Moreover determinism and laminarity (DET and 
LAM) are smaller smaller telling that that system is less regular. Consequently 
the more peculiar distribution of line lengths is confirmed by Lentr which is 
larger for the solution '2'. Additionally, in Fig. 10 we present the results of RR 
versus e (for relatively small e) . One can see the significant difference between both 
solutions. 

4 Summary and Conclusions 

We have analyzed the two DOF quarter-car model, assuming that damping, sus- 
pension through the unsprung mass excited by the road profile corrugation can 
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act as an perturbation on the main sprung mass. The obtained Melnikov criterion 
was latter confirmed by numerical simulations. The main conclusion coming from 
that point would be loss of stability of system appearing as the chaotic or tran- 
sient chaotic or escape solution. The present investigation is going beyond research 
dealing with a single DOF quarter-car model [4l ll8lll"9l l24]. In particular, the sin- 
gle DOF model assumes that the unsprung mass is significantly smaller than the 
sprung mass. 

One should note that the model used in this paper, however more realistic 
than previous single-degree-of-freedom ones, is relatively simple and would not be 
sufficient to simulate the detailed response of a vehicle or compare to experimental 
results from real vehicles. Unfortunately, more sophisticated half-car and full-car 
models [31,32,33 cannot be used in the frame of the presented approach as the 
heteroclinc trajectories could not be defined reliably. Furthermore, in higher DOF 
systems the analytic perturbation calculations are not possible. 

However, the present quarter-car model is able to capture the major nonlinear 
effects that occur in vehicle dynamics and has demonstrated the transition to 
chaotic vibrations and synchronization phenomena 12511321139] . Interestingly, the 
resulting critical amplitude curve (Fig. 4a, Eq. 19) has the maximum for w ss 3.1 
and the minimum for ui m 2.2. The minimum is obviously related to the resonance 
region of the decoupled unsprung solution X2 (Eqs. 8 and 9). 

It should be also noted that the recurrence plots technique appeared to be 
very useful to study the transient signals. Thus the conclusions came from the 
analytic approach have been confirmed. Indeed, this method is designed for the 
short time series [441149] Interestingly it also works for non-stationary signals [441 
|4"8] . The recurrences for RP and RQA analyzes have been obtain using the available 
commandline code written by Marwan [50] . The recurence approach can be also 
used to higher DOF models of vehicle dynamics. The corresponding results on half 
and full vehicle models will reported in a separate paper. 
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